Install / load packages needed:
data <- read_csv("data/cleanData.csv")
## Rows: 284 Columns: 131
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (9): id, sex, condition, Frage1, Frage2, Frage3, Leiter, Anmerkungen,...
## dbl (121): time, iat, ccs1, ccs2, ccs3, ccs4, ccs5, ccs6, ccs7, ccs8, ccs9,...
## dttm (1): StartDate
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
data <- data %>% dplyr::select(
id, time, iat, ccs, nr, nep, ipq, sod, ses, age, edu, sex, pol, vr_exp, vr_eval1, vr_eval2, vr_eval3,
vr_eval4, vr_eval5, span, seen, condition, starts_with("Frage"), hr_mean, Leiter, Anmerkungen, Zeit
)
head(data)
## # A tibble: 6 × 29
## id time iat ccs nr nep ipq sod ses age edu sex pol
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 2828… 1 0.179 1.08 3.48 4.07 NA NA 1.76 54 5 W 3
## 2 2828… 2 0.409 1 3.10 4.2 NA NA 1.71 54 5 W 3
## 3 7799… 1 -0.496 1.17 4.14 3.93 NA NA 1.53 21 4 M 2
## 4 7799… 2 -0.362 1.08 4.38 4.27 NA NA 1.35 21 4 M 2
## 5 4379… 1 0.517 1.33 3.43 3.73 NA NA 1.59 25 5 W 2
## 6 4379… 2 0.634 1.33 3.38 3.93 NA NA 1.53 25 5 W 2
## # … with 16 more variables: vr_exp <dbl>, vr_eval1 <dbl>, vr_eval2 <dbl>,
## # vr_eval3 <dbl>, vr_eval4 <dbl>, vr_eval5 <dbl>, span <dbl>, seen <dbl>,
## # condition <chr>, Frage1 <chr>, Frage2 <chr>, Frage3 <chr>, hr_mean <dbl>,
## # Leiter <chr>, Anmerkungen <chr>, Zeit <chr>
# factor for vr or not
data <- data %>% group_by(id) %>%
mutate(
vr = ifelse(condition %in% c("a", "b", "c"), TRUE, FALSE),
condition = factor(condition, levels = c("b", "a", "c", "video", "text.bild", "text"))
)
Keep in mind the conditions coding:
a == abstract
b == realistic
c == realistic but badly so
# mvoutlier::chisq.plot(data[,c(3:6)])
# [1] 124 162 161 29 54
# removing three most extreme cases
rmId <- data$id[c(124, 161, 162)]
data <- data %>% filter(!id %in% rmId)
data$id[c(142,144,275)]
## [1] "90811768" "78489627" "04901439"
We try to find an acceptable model for each DV.
First, I would like to calculate a principal component of all dependent variables (dvs)
df_env <- data[c("iat", "ccs", "nr", "nep")]
psych::fa.parallel(df_env)
## Parallel analysis suggests that the number of factors = 1 and the number of components = 1
prc_env <- princomp(df_env, cor = TRUE)
# Seems like one factor might just be enough. However, this may be more revealing when done
# on raw data
summary(prc_env)
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4
## Standard deviation 1.2980576 0.9835519 0.8829992 0.7536475
## Proportion of Variance 0.4212384 0.2418436 0.1949219 0.1419962
## Cumulative Proportion 0.4212384 0.6630819 0.8580038 1.0000000
data$env_pc <- prc_env$scores[,1]
Unidimensionality could be assumed. The scores of the first principal component were stored in data$env_pc. This vector can now be used as a dependent variable in further exploratory analyses.
This plot is not very useful I guess. Too crowded.
This section is concerned only with the VR conditions a, b, c. Specifically with the variables vr_eval 1:5 and with the sod and presence scale IPQ.
Check for outliers on these scales:
#check.data <- data %>% dplyr::filter((vr == T & time == 1 ) & !is.na(ipq))
check.data <- data %>% dplyr::filter((time == 1 ) & !is.na(ipq))
#outlier.data <- data %>% ungroup() %>% dplyr::filter(vr == T & time == 1) %>% dplyr::select(starts_with("vr_eval"), sod, ipq) %>%
# drop_na()
#mvoutlier::chisq.plot(check.data %>% ungroup() %>% dplyr::select(starts_with("vr_eval"), sod, ipq))
# remove: 38 1 63
# which corresponds to the ids:
#remove.ids <- check.data$id[c(38, 1, 63)]
# "44466757" "32504483" "80688810"
vars <-c("vr_eval1", "vr_eval2", "vr_eval3", "vr_eval4", "vr_eval5", "ipq", "sod")
desc_plot_data <- gather(data, specific, value, vars) %>%
# filter(!is.na(ipq)) %>%
arrange(id, specific) %>%
mutate(specific = ifelse(specific=="vr_eval1", "excitement",
ifelse(specific=="vr_eval2", "graphically pleasing",
ifelse(specific=="vr_eval3", "pleasant",
ifelse(specific=="vr_eval4", "realistic",
ifelse(specific=="vr_eval5", "enjoyment", specific))))),
VRE = ifelse(condition == "b", "R+",
ifelse(condition=="c", "R-",
ifelse(condition=="a", "A+", condition)))) %>%
dplyr::select(specific, value, id, VRE)
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(vars)` instead of `vars` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
(vr_eval_plot <- ggplot(data = desc_plot_data, aes(x = VRE, y = value, color = VRE)) +
facet_wrap( ~specific, nrow = 2)+
geom_violin(draw_quantiles = .5) + #, position = position_jitterdodge(dodge.width = 0.8, jitter.width = 0, jitter.height = 0))+
geom_point(alpha = .3, position = "jitter")+#, position = position_jitterdodge(dodge.width = 0, jitter.width = .2, jitter.height = .3)) +
ggthemes::theme_tufte() +
ylab("value") +
xlab("question")+
#labs(title="Evaluation of Virtual Environments")+
scale_y_continuous(breaks = c(1,3,5,7)) +
theme(legend.position = "bottom") )
## Warning: Removed 304 rows containing non-finite values (stat_ydensity).
## Warning: Removed 304 rows containing missing values (geom_point).
ipq_sod_plot_data <- gather(data, scale, value, c("ipq", "sod")) %>%
filter(time==1, condition %in% c("a", "b", "c", "video")) %>%
arrange(id, scale) %>%
mutate(VRE = ifelse(condition == "b", "R+",
ifelse(condition=="c", "R-",
ifelse(condition == "a", "A+", condition)))) %>%
dplyr::select(scale, value, id, VRE)
ipq_sod_plot <- ggplot(data = ipq_sod_plot_data, aes(x = VRE, y = value, color = VRE)) +
facet_wrap(~scale, ncol = 2)+
geom_violin(draw_quantiles = .5, position = position_jitterdodge(dodge.width = 0.8, jitter.width = 0, jitter.height = 0))+
geom_point(alpha = .3, position = position_jitterdodge(dodge.width = 0.8, jitter.width = .2, jitter.height = .1)) +
ggthemes::theme_tufte() +
ylab("scores") +
xlab("scale")+
#labs(title="Presence (IPQ) and Suspension of Disbelief (SOD)")+
scale_y_continuous(breaks = c(1,3,5,7), limits = c(1,7)); ipq_sod_plot
## Warning: Removed 4 rows containing non-finite values (stat_ydensity).
## Warning: Removed 6 rows containing missing values (geom_point).
So we will create two models for each dependent variables:
First, the model will have the formula:
dv ~ condition * time + (time | id)
This will be simplified to the following model if model fit is singular:
dv ~ condition * time + (1 | id)
This will estimate a random intercept for each participant. This model will only take as input the vr conditions (a, b & c), or: vr == TRUE.
The second model will have the formula:
dv ~ vr * time + (time | condition) + (1 | id)
Where a random slope for time is estimated per condition. Further, there is a random intercept per condition, and per id.
Should model fit be singular, we would simplify the model to:
dv ~ vr * time + (1 | condition) + (1 | id)
If still singular, we would simplify to:
dv ~ vr * time + (1 | id)
fit.lme <- function(form, dat){
lme4::lmer(formula = form, data = dat)
}
fit_models <- function(dv, dat){
# this function returns a function which fits a model based on a formula minus the predictors.
# This function can be used in the next function which implements the conditions for reducing model complexity if model fit is singular.
function(predictors){
form <- formula(paste(dv, predictors, sep = " ~ "))
print(form)
fit <- fit.lme(form = form, dat = dat)
}
}
predictors.vr <- c("condition * time + (time | id)", "condition * time + (1 | id)")
predictors.all <- c("vr * time + (time | condition) + (1 | id)", "vr * time + (-1 + time | condition) + (1 | id)", "vr * time + (1 | id)")
predictors.all2 <- c("vr * time + (time | condition) + (1 | id)", "vr * time + (1 | condition) + (1 | id)", "vr * time + (1 | id)")
# function to fit various models based on different inputs of predictors
fit_many <- function(pred.vector, dat, dv){
fit_model <- fit_models(dv, dat)
sing <- TRUE
i <- 1
while((sing) & i<=length(pred.vector)){
model <- try(fit_model(pred.vector[i]))
if(class(model)!="try-error"){
sing <- isSingular(model)
}
i <- i + 1
}
print(paste("is model singular: ", sing))
model
}
Vector containing name of all dv’s
dvs <- c("iat", "ccs", "nr", "nep", "env_pc")
# split data frame:
data.vr <- data %>% filter(vr)
vr.models <- lapply(dvs, FUN = function(dv) fit_many(pred.vector = predictors.vr, dat = data.vr, dv = dv))
## iat ~ condition * time + (time | id)
## <environment: 0x7fbb8206bbd8>
## Error : number of observations (=138) <= number of random effects (=138) for term (time | id); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
## iat ~ condition * time + (1 | id)
## <environment: 0x7fbb73d4aee8>
## [1] "is model singular: FALSE"
## ccs ~ condition * time + (time | id)
## <environment: 0x7fbb736b8180>
## Error : number of observations (=138) <= number of random effects (=138) for term (time | id); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
## ccs ~ condition * time + (1 | id)
## <environment: 0x7fbb7398ca88>
## [1] "is model singular: FALSE"
## nr ~ condition * time + (time | id)
## <environment: 0x7fbb8666b5e8>
## Error : number of observations (=138) <= number of random effects (=138) for term (time | id); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
## nr ~ condition * time + (1 | id)
## <environment: 0x7fbb71db7710>
## [1] "is model singular: FALSE"
## nep ~ condition * time + (time | id)
## <environment: 0x7fbb70714d20>
## Error : number of observations (=138) <= number of random effects (=138) for term (time | id); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
## nep ~ condition * time + (1 | id)
## <environment: 0x7fbb7188a548>
## [1] "is model singular: FALSE"
## env_pc ~ condition * time + (time | id)
## <environment: 0x7fbb70905f00>
## Error : number of observations (=138) <= number of random effects (=138) for term (time | id); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
## env_pc ~ condition * time + (1 | id)
## <environment: 0x7fbb9067f070>
## [1] "is model singular: FALSE"
all.models <- lapply(dvs, FUN = function(dv) fit_many(pred.vector = predictors.all, dat = data, dv = dv))
## iat ~ vr * time + (time | condition) + (1 | id)
## <environment: 0x7fbb87c40500>
## boundary (singular) fit: see ?isSingular
## iat ~ vr * time + (-1 + time | condition) + (1 | id)
## <environment: 0x7fbb83f83b08>
## [1] "is model singular: FALSE"
## ccs ~ vr * time + (time | condition) + (1 | id)
## <environment: 0x7fbb85df3878>
## boundary (singular) fit: see ?isSingular
## ccs ~ vr * time + (-1 + time | condition) + (1 | id)
## <environment: 0x7fbb83691bf8>
## boundary (singular) fit: see ?isSingular
## ccs ~ vr * time + (1 | id)
## <environment: 0x7fbb84ae3668>
## [1] "is model singular: FALSE"
## nr ~ vr * time + (time | condition) + (1 | id)
## <environment: 0x7fbb859b00d0>
## boundary (singular) fit: see ?isSingular
## nr ~ vr * time + (-1 + time | condition) + (1 | id)
## <environment: 0x7fbba3f73380>
## [1] "is model singular: FALSE"
## nep ~ vr * time + (time | condition) + (1 | id)
## <environment: 0x7fbb72d85c80>
## boundary (singular) fit: see ?isSingular
## nep ~ vr * time + (-1 + time | condition) + (1 | id)
## <environment: 0x7fbb752992d0>
## [1] "is model singular: FALSE"
## env_pc ~ vr * time + (time | condition) + (1 | id)
## <environment: 0x7fbb7370f910>
## boundary (singular) fit: see ?isSingular
## env_pc ~ vr * time + (-1 + time | condition) + (1 | id)
## <environment: 0x7fbb70bc3c48>
## boundary (singular) fit: see ?isSingular
## env_pc ~ vr * time + (1 | id)
## <environment: 0x7fbb71f874e0>
## [1] "is model singular: FALSE"
I save model diagnostics as pdfs separately, for visibility reasons.
plot_diagn <- function(model){
filename <- paste( model@call$formula[2], sub("\\ .*", "", model@call$formula[3]), sep = "_")
png(filename = paste("analysisOutputs/diagnostics/", filename, ".png", sep = ""), # The directory you want to save the file in
#paper = "a3",
height = 5900/4,
width = 4200/4
)
print(performance::check_model(model) )
dev.off()
}
lapply(vr.models, FUN = plot_diagn)
## [[1]]
## quartz_off_screen
## 2
##
## [[2]]
## quartz_off_screen
## 2
##
## [[3]]
## quartz_off_screen
## 2
##
## [[4]]
## quartz_off_screen
## 2
##
## [[5]]
## quartz_off_screen
## 2
lapply(all.models, FUN = plot_diagn)
## [[1]]
## quartz_off_screen
## 2
##
## [[2]]
## quartz_off_screen
## 2
##
## [[3]]
## quartz_off_screen
## 2
##
## [[4]]
## quartz_off_screen
## 2
##
## [[5]]
## quartz_off_screen
## 2
I focus model diagnostic on the vr models. They include all data. Residuals are slightly left skewed. However, this does not yet warrant a transformation of the dv in my opinion.
iat_vr_diagnostics
Some thoughts: Band of residuals increases as fitted values increase. Homogeneity of variance seems acceptable. Random effects appear normal.
ccs_vr_diagnostics
Some thoughts: Homogeneity of variance appears implausible. Residual variance increases with larger fitted values.
Residuals are also not normally distributed. Random effects do not appear normal.
The reason for this unexpected behaviour may well be the floor-effect of the dv ccs. There was generally a very low ccs score for participants. This is due to the relatively extreme nature of climate change scepticism, especially in a relatively well educated sample.
Maybe a boxcox transformation may help:
#estimate lambda of the boxcox transformation
bc <- boxcox(ccs ~ vr * time, data = data)
lambda_ccs <- bc$x[which.max(bc$y)]
# transform data according to the transformation
data <- data %>%
mutate(ccs_bc = (ccs^lambda_ccs-1)/lambda_ccs)
# refit the model
all.ccs2 <- fit_many(pred.vector = predictors.all, dat = data, dv = "ccs_bc")
## ccs_bc ~ vr * time + (time | condition) + (1 | id)
## <environment: 0x7fbb729cdd58>
## boundary (singular) fit: see ?isSingular
## ccs_bc ~ vr * time + (-1 + time | condition) + (1 | id)
## <environment: 0x7fbb75293078>
## boundary (singular) fit: see ?isSingular
## ccs_bc ~ vr * time + (1 | id)
## <environment: 0x7fbb70ad7618>
## [1] "is model singular: FALSE"
performance::check_model(all.ccs2)
The situation has improved! All model assumptions appear plausible.
all.models[[2]] <- all.ccs2
For within the VE:
ccs_condition_diagnostics
And based on the transformed ccs:
data.vr <- data.vr %>%
mutate(ccs_bc = (ccs^lambda_ccs-1)/lambda_ccs)
vr.ccs2 <- fit_many(pred.vector = predictors.vr, dat = data.vr, dv = "ccs_bc")
## ccs_bc ~ condition * time + (time | id)
## <environment: 0x7fbb57caa4b8>
## Error : number of observations (=138) <= number of random effects (=138) for term (time | id); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable
## ccs_bc ~ condition * time + (1 | id)
## <environment: 0x7fbb70126b88>
## [1] "is model singular: FALSE"
performance::check_model(vr.ccs2)
vr.models[[2]] <- vr.ccs2
nr_vr_diagnostics
Model assumptions are not too far off: Slight slope in the “fitted vs residuals”. Homogeneity assumption is appropriate. Residual distribution is slightly skewed with a heavy left tail. ID intercept distribution is not quite normal, but not far from it. Slightly skewed as well.
Overall assumptions seem acceptable and warrant no further action.
nep_vr_diagnostics
Model assumptions are not too far off: Slight slope in the “fitted vs residuals”. Residuals appear normally distributed, slightly skewed to the right.
Overall assumptions seem acceptable and warrant no further action.
env_pc_vr_diagnostics
min(data$env_pc)
## [1] -4.954711
data$env_pc2 <- data$env_pc + 6
#estimate lambda of the boxcox transformation
bc <- boxcox(env_pc2 ~ vr * time, data = data)
lambda_pc <- bc$x[which.max(bc$y)]
# transform data according to the transformation
data <- data %>%
mutate(env_pc_bc = (env_pc2^lambda_pc-1)/lambda_pc)
# refit the model
all.env_pc2 <- fit_many(pred.vector = predictors.all, dat = data, dv = "env_pc_bc")
## env_pc_bc ~ vr * time + (time | condition) + (1 | id)
## <environment: 0x7fbb900d0758>
## boundary (singular) fit: see ?isSingular
## env_pc_bc ~ vr * time + (-1 + time | condition) + (1 | id)
## <environment: 0x7fbb7435a608>
## boundary (singular) fit: see ?isSingular
## env_pc_bc ~ vr * time + (1 | id)
## <environment: 0x7fbb7078ca80>
## [1] "is model singular: FALSE"
performance::check_model(all.env_pc2)
Transformation does not help much here. Lets not do it.
But let us save the final model diagnostics plots as well.
lapply(vr.models, FUN = plot_diagn)
## [[1]]
## quartz_off_screen
## 2
##
## [[2]]
## quartz_off_screen
## 2
##
## [[3]]
## quartz_off_screen
## 2
##
## [[4]]
## quartz_off_screen
## 2
##
## [[5]]
## quartz_off_screen
## 2
lapply(all.models, FUN = plot_diagn)
## [[1]]
## quartz_off_screen
## 2
##
## [[2]]
## quartz_off_screen
## 2
##
## [[3]]
## quartz_off_screen
## 2
##
## [[4]]
## quartz_off_screen
## 2
##
## [[5]]
## quartz_off_screen
## 2
lapply(vr.models, FUN = lmerTest:::summary.lmerModLmerTest)
## Coercing object to class 'lmerModLmerTest'
## Coercing object to class 'lmerModLmerTest'
## Coercing object to class 'lmerModLmerTest'
## Coercing object to class 'lmerModLmerTest'
## Coercing object to class 'lmerModLmerTest'
## [[1]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: iat ~ condition * time + (1 | id)
## Data: dat
##
## REML criterion at convergence: 173.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.60384 -0.59622 -0.01456 0.55523 1.93986
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.1132 0.3365
## Residual 0.1068 0.3268
## Number of obs: 138, groups: id, 69
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.17718 0.16773 106.85910 1.056 0.293
## conditiona -0.32465 0.23721 106.85910 -1.369 0.174
## conditionc 0.02806 0.23721 106.85910 0.118 0.906
## time -0.01191 0.09636 66.00000 -0.124 0.902
## conditiona:time 0.11634 0.13627 66.00000 0.854 0.396
## conditionc:time 0.11535 0.13627 66.00000 0.846 0.400
##
## Correlation of Fixed Effects:
## (Intr) condtn cndtnc time cndtn:t
## conditiona -0.707
## conditionc -0.707 0.500
## time -0.862 0.609 0.609
## conditin:tm 0.609 -0.862 -0.431 -0.707
## conditnc:tm 0.609 -0.431 -0.862 -0.707 0.500
##
## [[2]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ccs_bc ~ condition * time + (1 | id)
## Data: dat
##
## REML criterion at convergence: -153.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.88848 -0.39218 -0.03104 0.48900 1.99181
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.022653 0.15051
## Residual 0.004976 0.07054
## Number of obs: 138, groups: id, 69
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.747e-01 4.546e-02 1.316e+02 3.844 0.000188 ***
## conditiona 7.150e-02 6.429e-02 1.316e+02 1.112 0.268108
## conditionc 1.003e-01 6.429e-02 1.316e+02 1.561 0.120977
## time 4.354e-04 2.080e-02 6.600e+01 0.021 0.983363
## conditiona:time -8.028e-03 2.942e-02 6.600e+01 -0.273 0.785799
## conditionc:time -3.700e-02 2.942e-02 6.600e+01 -1.258 0.212892
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) condtn cndtnc time cndtn:t
## conditiona -0.707
## conditionc -0.707 0.500
## time -0.686 0.485 0.485
## conditin:tm 0.485 -0.686 -0.343 -0.707
## conditnc:tm 0.485 -0.343 -0.686 -0.707 0.500
##
## [[3]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: nr ~ condition * time + (1 | id)
## Data: dat
##
## REML criterion at convergence: 98.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.74225 -0.50265 -0.03626 0.54439 1.71416
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.11846 0.3442
## Residual 0.04125 0.2031
## Number of obs: 138, groups: id, 69
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.90269 0.11882 129.34427 32.847 <2e-16 ***
## conditiona -0.38509 0.16803 129.34427 -2.292 0.0235 *
## conditionc 0.09731 0.16803 129.34427 0.579 0.5635
## time 0.09110 0.05989 66.00000 1.521 0.1330
## conditiona:time 0.02277 0.08470 66.00000 0.269 0.7888
## conditionc:time -0.06625 0.08470 66.00000 -0.782 0.4369
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) condtn cndtnc time cndtn:t
## conditiona -0.707
## conditionc -0.707 0.500
## time -0.756 0.535 0.535
## conditin:tm 0.535 -0.756 -0.378 -0.707
## conditnc:tm 0.535 -0.378 -0.756 -0.707 0.500
##
## [[4]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: nep ~ condition * time + (1 | id)
## Data: dat
##
## REML criterion at convergence: 93.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.32935 -0.59590 -0.04369 0.49713 1.85313
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.13842 0.3721
## Residual 0.03437 0.1854
## Number of obs: 138, groups: id, 69
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.675362 0.116144 131.998743 31.645 <2e-16 ***
## conditiona -0.014493 0.164253 131.998743 -0.088 0.9298
## conditionc 0.150725 0.164253 131.998743 0.918 0.3605
## time 0.115942 0.054666 66.000003 2.121 0.0377 *
## conditiona:time 0.008696 0.077310 66.000003 0.112 0.9108
## conditionc:time -0.078261 0.077310 66.000003 -1.012 0.3151
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) condtn cndtnc time cndtn:t
## conditiona -0.707
## conditionc -0.707 0.500
## time -0.706 0.499 0.499
## conditin:tm 0.499 -0.706 -0.353 -0.707
## conditnc:tm 0.499 -0.353 -0.706 -0.707 0.500
##
## [[5]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: env_pc ~ condition * time + (1 | id)
## Data: dat
##
## REML criterion at convergence: 367.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61435 -0.42168 0.01546 0.43108 1.86555
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.2457 1.1161
## Residual 0.2474 0.4974
## Number of obs: 138, groups: id, 69
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.16537 0.32856 130.61147 -0.503 0.6156
## conditiona -0.74721 0.46465 130.61147 -1.608 0.1102
## conditionc 0.10299 0.46465 130.61147 0.222 0.8249
## time 0.25199 0.14668 66.00000 1.718 0.0905 .
## conditiona:time 0.07992 0.20744 66.00000 0.385 0.7013
## conditionc:time -0.05075 0.20744 66.00000 -0.245 0.8075
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) condtn cndtnc time cndtn:t
## conditiona -0.707
## conditionc -0.707 0.500
## time -0.670 0.474 0.474
## conditin:tm 0.474 -0.670 -0.335 -0.707
## conditnc:tm 0.474 -0.335 -0.670 -0.707 0.500
contrast analysis:
# iat
contrast( emmeans(vr.models[[1]], ~ time | condition))
## condition = b:
## contrast estimate SE df t.ratio p.value
## 1 effect 0.00595 0.0482 66 0.124 0.9020
## 2 effect -0.00595 0.0482 66 -0.124 0.9020
##
## condition = a:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.05221 0.0482 66 -1.084 0.2824
## 2 effect 0.05221 0.0482 66 1.084 0.2824
##
## condition = c:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.05172 0.0482 66 -1.073 0.2870
## 2 effect 0.05172 0.0482 66 1.073 0.2870
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: fdr method for 2 tests
# ccs
contrast( emmeans(vr.models[[2]], ~ time | condition))
## condition = b:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.000218 0.0104 66 -0.021 0.9834
## 2 effect 0.000218 0.0104 66 0.021 0.9834
##
## condition = a:
## contrast estimate SE df t.ratio p.value
## 1 effect 0.003796 0.0104 66 0.365 0.7163
## 2 effect -0.003796 0.0104 66 -0.365 0.7163
##
## condition = c:
## contrast estimate SE df t.ratio p.value
## 1 effect 0.018284 0.0104 66 1.758 0.0834
## 2 effect -0.018284 0.0104 66 -1.758 0.0834
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: fdr method for 2 tests
# nr
contrast( emmeans(vr.models[[3]], ~ time | condition))
## condition = b:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.0455 0.0299 66 -1.521 0.1330
## 2 effect 0.0455 0.0299 66 1.521 0.1330
##
## condition = a:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.0569 0.0299 66 -1.901 0.0616
## 2 effect 0.0569 0.0299 66 1.901 0.0616
##
## condition = c:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.0124 0.0299 66 -0.415 0.6796
## 2 effect 0.0124 0.0299 66 0.415 0.6796
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: fdr method for 2 tests
# nep
contrast( emmeans(vr.models[[4]], ~ time | condition))
## condition = b:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.0580 0.0273 66 -2.121 0.0377
## 2 effect 0.0580 0.0273 66 2.121 0.0377
##
## condition = a:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.0623 0.0273 66 -2.280 0.0258
## 2 effect 0.0623 0.0273 66 2.280 0.0258
##
## condition = c:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.0188 0.0273 66 -0.689 0.4931
## 2 effect 0.0188 0.0273 66 0.689 0.4931
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: fdr method for 2 tests
lapply(all.models, FUN = lmerTest:::summary.lmerModLmerTest)
## Coercing object to class 'lmerModLmerTest'
## Coercing object to class 'lmerModLmerTest'
## Coercing object to class 'lmerModLmerTest'
## Coercing object to class 'lmerModLmerTest'
## Coercing object to class 'lmerModLmerTest'
## [[1]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: iat ~ vr * time + (-1 + time | condition) + (1 | id)
## Data: dat
##
## REML criterion at convergence: 312.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.00917 -0.60925 0.02082 0.59421 2.05824
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.098991 0.31463
## condition time 0.002013 0.04486
## Residual 0.094707 0.30775
## Number of obs: 284, groups: id, 142; condition, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.38251 0.08856 225.15572 4.319 2.35e-05 ***
## vrTRUE -0.30420 0.12704 225.15572 -2.394 0.0175 *
## time -0.08398 0.05715 19.60552 -1.470 0.1575
## vrTRUE:time 0.14930 0.08174 20.39708 1.827 0.0824 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vrTRUE time
## vrTRUE -0.697
## time -0.769 0.536
## vrTRUE:time 0.538 -0.771 -0.699
##
## [[2]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ccs_bc ~ vr * time + (1 | id)
## Data: dat
##
## REML criterion at convergence: -351.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.33016 -0.44377 -0.04596 0.56783 2.12335
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.018962 0.13770
## Residual 0.005645 0.07513
## Number of obs: 284, groups: id, 142
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.236218 0.025424 278.363625 9.291 <2e-16 ***
## vrTRUE -0.004193 0.036472 278.363625 -0.115 0.909
## time -0.008953 0.012436 140.000002 -0.720 0.473
## vrTRUE:time -0.005622 0.017840 140.000002 -0.315 0.753
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vrTRUE time
## vrTRUE -0.697
## time -0.734 0.511
## vrTRUE:time 0.511 -0.734 -0.697
##
## [[3]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: nr ~ vr * time + (-1 + time | condition) + (1 | id)
## Data: dat
##
## REML criterion at convergence: 262.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.59554 -0.44842 0.02171 0.48568 1.76425
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.2263312 0.47574
## condition time 0.0003591 0.01895
## Residual 0.0400008 0.20000
## Number of obs: 284, groups: id, 142; condition, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.828441 0.076421 268.099977 50.097 <2e-16 ***
## vrTRUE -0.021678 0.109631 268.099978 -0.198 0.843
## time -0.001225 0.034866 9.120707 -0.035 0.973
## vrTRUE:time 0.077830 0.049948 9.573006 1.558 0.152
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vrTRUE time
## vrTRUE -0.697
## time -0.617 0.430
## vrTRUE:time 0.431 -0.618 -0.698
##
## [[4]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: nep ~ vr * time + (-1 + time | condition) + (1 | id)
## Data: dat
##
## REML criterion at convergence: 184.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.13567 -0.52383 -0.01767 0.49178 2.18200
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.1427876 0.37787
## condition time 0.0001482 0.01217
## Residual 0.0352898 0.18786
## Number of obs: 284, groups: id, 142; condition, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.79087 0.06613 276.21351 57.325 <2e-16 ***
## vrTRUE -0.07009 0.09487 276.21351 -0.739 0.461
## time 0.05487 0.03188 14.13761 1.721 0.107
## vrTRUE:time 0.03788 0.04570 14.88848 0.829 0.420
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vrTRUE time
## vrTRUE -0.697
## time -0.688 0.480
## vrTRUE:time 0.480 -0.688 -0.698
##
## [[5]]
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: env_pc ~ vr * time + (1 | id)
## Data: dat
##
## REML criterion at convergence: 766.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.33852 -0.41115 0.01461 0.44995 1.90160
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.4722 1.2133
## Residual 0.2278 0.4773
## Number of obs: 284, groups: id, 142
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.09543 0.18914 267.65387 -0.505 0.6143
## vrTRUE -0.28468 0.27133 267.65387 -1.049 0.2950
## time 0.05577 0.07901 140.00000 0.706 0.4815
## vrTRUE:time 0.20595 0.11334 140.00000 1.817 0.0713 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) vrTRUE time
## vrTRUE -0.697
## time -0.627 0.437
## vrTRUE:time 0.437 -0.627 -0.697
Follow up tests:
# iat
contrast( emmeans(all.models[[1]], ~ time | vr))
## vr = FALSE:
## contrast estimate SE df t.ratio p.value
## 1 effect 0.0420 0.0286 21.2 1.470 0.1563
## 2 effect -0.0420 0.0286 21.2 -1.470 0.1563
##
## vr = TRUE:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.0327 0.0292 23.0 -1.118 0.2753
## 2 effect 0.0327 0.0292 23.0 1.118 0.2753
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: fdr method for 2 tests
# ccs
contrast( emmeans(all.models[[2]], ~ time | vr))
## vr = FALSE:
## contrast estimate SE df t.ratio p.value
## 1 effect 0.00448 0.00622 140 0.720 0.4728
## 2 effect -0.00448 0.00622 140 -0.720 0.4728
##
## vr = TRUE:
## contrast estimate SE df t.ratio p.value
## 1 effect 0.00729 0.00640 140 1.139 0.2565
## 2 effect -0.00729 0.00640 140 -1.139 0.2565
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: fdr method for 2 tests
# nr
contrast( emmeans(all.models[[3]], ~ time | vr))
## vr = FALSE:
## contrast estimate SE df t.ratio p.value
## 1 effect 0.000613 0.0174 9.67 0.035 0.9727
## 2 effect -0.000613 0.0174 9.67 -0.035 0.9727
##
## vr = TRUE:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.038302 0.0179 10.66 -2.142 0.0562
## 2 effect 0.038302 0.0179 10.66 2.142 0.0562
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: fdr method for 2 tests
# nep
contrast( emmeans(all.models[[4]], ~ time | vr))
## vr = FALSE:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.0274 0.0159 13.0 -1.721 0.1090
## 2 effect 0.0274 0.0159 13.0 1.721 0.1090
##
## vr = TRUE:
## contrast estimate SE df t.ratio p.value
## 1 effect -0.0464 0.0164 14.4 -2.833 0.0130
## 2 effect 0.0464 0.0164 14.4 2.833 0.0130
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: fdr method for 2 tests
For the three VE conditions we need ot test the predictor condition with model comparisons.
compare.models <- function(model){
model0 <- update(model, .~. - time:condition)
anova(model0, model)
}
lapply(vr.models, compare.models)
## refitting model(s) with ML (instead of REML)
## refitting model(s) with ML (instead of REML)
## refitting model(s) with ML (instead of REML)
## refitting model(s) with ML (instead of REML)
## refitting model(s) with ML (instead of REML)
## [[1]]
## Data: dat
## Models:
## model0: iat ~ condition + time + (1 | id)
## model: iat ~ condition * time + (1 | id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model0 6 168.31 185.87 -78.155 156.31
## model 8 171.31 194.73 -77.655 155.31 1.0001 2 0.6065
##
## [[2]]
## Data: dat
## Models:
## model0: ccs_bc ~ condition + time + (1 | id)
## model: ccs_bc ~ condition * time + (1 | id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model0 6 -172.93 -155.37 92.464 -184.93
## model 8 -170.74 -147.32 93.368 -186.74 1.8068 2 0.4052
##
## [[3]]
## Data: dat
## Models:
## model0: nr ~ condition + time + (1 | id)
## model: nr ~ condition * time + (1 | id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model0 6 90.456 108.02 -39.228 78.456
## model 8 93.220 116.64 -38.610 77.220 1.2358 2 0.5391
##
## [[4]]
## Data: dat
## Models:
## model0: nep ~ condition + time + (1 | id)
## model: nep ~ condition * time + (1 | id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model0 6 85.962 103.53 -36.981 73.962
## model 8 88.376 111.79 -36.188 72.376 1.5864 2 0.4524
##
## [[5]]
## Data: dat
## Models:
## model0: env_pc ~ condition + time + (1 | id)
## model: env_pc ~ condition * time + (1 | id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model0 6 371.07 388.63 -179.53 359.07
## model 8 374.64 398.06 -179.32 358.64 0.4205 2 0.8104
The condition:time interaction did not significantly add to the explained variance.
First, we define a couple of helping functions
make.x.table <- function(lmermods, save = TRUE, add = ""){
lapply(lmermods, function(x) formula(x)[2])
dvs <- as.character(lapply(lmermods, function(x) as.character(formula(x)[2])))
iv <- sub("\\ .*", "", lmermods[[1]]@call$formula[3])
helpf <- function(x){
# get fixed effects table
fix.temp <- round(as.data.frame(summary(x)$coefficients), 3)
# calculate and round confidence intervals
conf.temp <- as.data.frame(confint(x, method = "profile"))
conf.temp <- round(conf.temp, 3)
conf.temp$CI <- paste("[", conf.temp[,1], "; ", conf.temp[,2],"]", sep = "")
# put together
fix.temp$CI <- conf.temp$CI[!startsWith(rownames(conf.temp), ".")]
fix.temp["DV"] <- ""
#fix.temp[paste(dv, " ~", sep = "")][1,] <-paste(formula(x))[3]
fix.temp["DV"][1,] <- "deleteme"
fix.temp <- rownames_to_column(fix.temp, var = "Coef")
fix.temp <- fix.temp[, c(6, 1, 2, 5, 3, 4)]
# return value
return(fix.temp)
}
# confidence intervals:
fixed.effects.list <- lapply(lmermods, FUN = helpf)
fixed.effects.df <- do.call("rbind", fixed.effects.list)
fixed.effects.df$DV[fixed.effects.df$DV=="deleteme"] <- dvs
#l.temp <- length(rownames(fixed.effects.df))
#cond <- toupper(strsplit(rownames(fixed.effects.df)[l.temp], split = " ")[[1]][5])
rownames(fixed.effects.df) <- NULL
tabl <- xtable(fixed.effects.df[-6],
caption = paste(toupper(iv), "Models", add),
label = paste("tab:",add , iv, "-models", sep = ""))
if(save){
print.xtable(tabl, file = paste("analysisOutputs/tables/", add, iv, "-modeltable.tex", sep = ""),
include.rownames=FALSE,
hline.after = c(-1, c(which(fixed.effects.df[1]!="")-1), nrow(tabl))#,
#add.to.row = list(list(-1), c(paste("\\hspace{-3mm}", toupper(dv), "$\\sim$%")))
)
}else{
print.xtable(tabl, #file = paste("analysisOutputs/tables/", add, dv, "-by", cond, "-modeltable.tex", sep = ""),
include.rownames=FALSE,
hline.after = c(-1, c(which(fixed.effects.df[1]!="")-1), nrow(tabl))#,
#add.to.row = list(list(-1), c(paste("\\hspace{-3mm}", toupper(dv), "$\\sim$%")))
)}
# which(fixed.effects.df[2]!="")-1
# this saves directly into my .tex folder
}
Then we create and save the tables:
make.x.table(all.models, save = FALSE)
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): Last two rows have
## identical or NA .zeta values: using minstep
## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): Last two rows have
## identical or NA .zeta values: using minstep
## Warning in zetafun(np, ns): slightly lower deviances (diff=-2.27374e-13)
## detected
## Warning in FUN(X[[i]], ...): non-monotonic profile for .sig02
## Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit
## for .sig02: falling back to linear interpolation
## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
## Computing profile confidence intervals ...
## % latex table generated in R 4.1.2 by xtable 1.8-4 package
## % Sat Feb 12 00:37:16 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{llrlr}
## \hline
## DV & Coef & Estimate & CI & Std. Error \\
## \hline
## iat & (Intercept) & 0.38 & [0.209; 0.556] & 0.09 \\
## & vrTRUE & -0.30 & [-0.553; -0.056] & 0.13 \\
## & time & -0.08 & [-0.19; 0.022] & 0.06 \\
## & vrTRUE:time & 0.15 & [-0.002; 0.301] & 0.08 \\
## \hline
## ccs\_bc & (Intercept) & 0.24 & [0.187; 0.286] & 0.03 \\
## & vrTRUE & -0.00 & [-0.075; 0.067] & 0.04 \\
## & time & -0.01 & [-0.033; 0.015] & 0.01 \\
## & vrTRUE:time & -0.01 & [-0.041; 0.029] & 0.02 \\
## \hline
## nr & (Intercept) & 3.83 & [3.679; 3.978] & 0.08 \\
## & vrTRUE & -0.02 & [-0.236; 0.193] & 0.11 \\
## & time & -0.00 & [-0.067; 0.065] & 0.04 \\
## & vrTRUE:time & 0.08 & [-0.017; 0.172] & 0.05 \\
## \hline
## nep & (Intercept) & 3.79 & [3.662; 3.92] & 0.07 \\
## & vrTRUE & -0.07 & [-0.255; 0.115] & 0.10 \\
## & time & 0.06 & [-0.006; 0.116] & 0.03 \\
## & vrTRUE:time & 0.04 & [-0.049; 0.125] & 0.05 \\
## \hline
## env\_pc & (Intercept) & -0.10 & [-0.465; 0.274] & 0.19 \\
## & vrTRUE & -0.28 & [-0.815; 0.245] & 0.27 \\
## & time & 0.06 & [-0.099; 0.211] & 0.08 \\
## & vrTRUE:time & 0.21 & [-0.016; 0.428] & 0.11 \\
## \hline
## \end{tabular}
## \caption{VR Models }
## \label{tab:vr-models}
## \end{table}
make.x.table(vr.models, save = FALSE)
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## % latex table generated in R 4.1.2 by xtable 1.8-4 package
## % Sat Feb 12 00:37:19 2022
## \begin{table}[ht]
## \centering
## \begin{tabular}{llrlr}
## \hline
## DV & Coef & Estimate & CI & Std. Error \\
## \hline
## iat & (Intercept) & 0.18 & [-0.147; 0.502] & 0.17 \\
## & conditiona & -0.33 & [-0.783; 0.134] & 0.24 \\
## & conditionc & 0.03 & [-0.431; 0.487] & 0.24 \\
## & time & -0.01 & [-0.199; 0.175] & 0.10 \\
## & conditiona:time & 0.12 & [-0.149; 0.381] & 0.14 \\
## & conditionc:time & 0.12 & [-0.15; 0.38] & 0.14 \\
## \hline
## ccs\_bc & (Intercept) & 0.17 & [0.087; 0.262] & 0.04 \\
## & conditiona & 0.07 & [-0.053; 0.196] & 0.06 \\
## & conditionc & 0.10 & [-0.024; 0.224] & 0.06 \\
## & time & 0.00 & [-0.04; 0.041] & 0.02 \\
## & conditiona:time & -0.01 & [-0.065; 0.049] & 0.03 \\
## & conditionc:time & -0.04 & [-0.094; 0.02] & 0.03 \\
## \hline
## nr & (Intercept) & 3.90 & [3.673; 4.132] & 0.12 \\
## & conditiona & -0.39 & [-0.709; -0.061] & 0.17 \\
## & conditionc & 0.10 & [-0.227; 0.422] & 0.17 \\
## & time & 0.09 & [-0.025; 0.208] & 0.06 \\
## & conditiona:time & 0.02 & [-0.142; 0.187] & 0.09 \\
## & conditionc:time & -0.07 & [-0.231; 0.098] & 0.09 \\
## \hline
## nep & (Intercept) & 3.67 & [3.451; 3.9] & 0.12 \\
## & conditiona & -0.01 & [-0.332; 0.303] & 0.16 \\
## & conditionc & 0.15 & [-0.166; 0.468] & 0.16 \\
## & time & 0.12 & [0.01; 0.222] & 0.06 \\
## & conditiona:time & 0.01 & [-0.142; 0.159] & 0.08 \\
## & conditionc:time & -0.08 & [-0.229; 0.072] & 0.08 \\
## \hline
## env\_pc & (Intercept) & -0.17 & [-0.8; 0.469] & 0.33 \\
## & conditiona & -0.75 & [-1.644; 0.15] & 0.47 \\
## & conditionc & 0.10 & [-0.794; 1] & 0.47 \\
## & time & 0.25 & [-0.033; 0.537] & 0.15 \\
## & conditiona:time & 0.08 & [-0.323; 0.483] & 0.21 \\
## & conditionc:time & -0.05 & [-0.454; 0.352] & 0.21 \\
## \hline
## \end{tabular}
## \caption{CONDITION Models }
## \label{tab:condition-models}
## \end{table}
make.x.table(all.models, save = TRUE)
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): Last two rows have
## identical or NA .zeta values: using minstep
## Warning in nextpar(mat, cc, i, delta, lowcut, upcut): Last two rows have
## identical or NA .zeta values: using minstep
## Warning in zetafun(np, ns): slightly lower deviances (diff=-2.27374e-13)
## detected
## Warning in FUN(X[[i]], ...): non-monotonic profile for .sig02
## Warning in confint.thpr(pp, level = level, zeta = zeta): bad spline fit
## for .sig02: falling back to linear interpolation
## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values
## Computing profile confidence intervals ...
make.x.table(vr.models, save = TRUE)
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
## Computing profile confidence intervals ...
I would like to get a table / data structure that gives me the mean values for each combination of relevant time and condition ### CI plot helping functions
ci.plot <- function(colorby = "condition"){
#colorby = "condition"
#df <- vr.res.md.comp.mean
pos <- ifelse(colorby=="condition", "dodge2", "identity")
function(df){
df$cond <- df[[colorby]]
if(colorby != "condition"){
df <- df %>% dplyr::filter(condition=="a"|condition=="text" )
}
ggplot(df, aes(x = time, y = ml.value, color = cond)) +
facet_wrap(~dv, scales = "free") +
geom_line(position = position_jitterdodge(dodge.width = 0.2, jitter.width = 0, jitter.height = 0)) +
geom_point(position = position_jitterdodge(dodge.width = 0.2, jitter.width = 0, jitter.height = 0)) +
geom_errorbar(aes(ymin = df$`2.5 %`, ymax = df$`97.5 %`), position = "dodge2", width = 0.2) +
ggthemes::theme_tufte() +
scale_x_discrete() +
xlim(c("before", "after"))
}
}
ci.plot.vr <- ci.plot("vr")
ci.plot.condition <- ci.plot("condition")
#ci.plot.condition(vr.res.md.comp.mean)
#ci.plot.vr(vr.nonvr.mean)
cond.mean.ci <- function(data.temp){
predict.fun <- function(my.lmm) {
# my.lmm <- md
#data.temp <- df.predicted.vr
predict(my.lmm, newdata = data.temp, re.form = NA) # This is predict.merMod
#data.temp$x <- x
#data.temp$y <- rep(x[data.temp$time==1]-x[data.temp$time==2], each = 2)
#y <- x
}
function(md){
data.temp$ml.value <- predict.fun(md)
# Make predictions in 100 bootstraps of the LMM. Use these to get confidence
# intervals.
dv <- as.character(formula(md))[2]
lmm.boots <- bootMer(md, predict.fun, nsim = 1000,
#parallel = "multicore",
ncpus = (detectCores(all.tests = FALSE, logical = TRUE)-1),
type = "semiparametric",
use.u = T)
data.temp <- cbind("dv" = dv, data.temp, confint(lmm.boots))
return(data.temp)
}
}
df.predicted.all <- data.frame(
id = rep(as.factor(1:6), each = 2),
time = rep(1:2, times = 6),
condition = rep(levels(as.factor(data$condition)), each = 2),
vr = rep(c(T, F), each = 6)
)
df.predicted.vr <- df.predicted.all %>% filter(vr)
Create function for each option:
cond.mean.ci.all <- cond.mean.ci(data = df.predicted.all)
cond.mean.ci.vr <- cond.mean.ci(data = df.predicted.vr)
First comparing VE to control conditions
# Bootstrap CI's
vr.nonvr.mean.list <- lapply(all.models[1:4], cond.mean.ci.all)
vr.nonvr.mean <- do.call("rbind", vr.nonvr.mean.list)
(vr.ci.plot <- ci.plot.vr(vr.nonvr.mean) +
guides(color=guide_legend(title="VR")))
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Use of `df$`2.5 %`` is discouraged. Use `2.5 %` instead.
## Warning: Use of `df$`97.5 %`` is discouraged. Use `97.5 %` instead.
Then for only between the VE conditions:
onlyvr.cond.mean.list <- lapply(vr.models[1:4], cond.mean.ci.vr)
onlyvr.cond.mean <- do.call("rbind", onlyvr.cond.mean.list)
(cond.ci.plot <- ci.plot.condition(onlyvr.cond.mean))
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Use of `df$`2.5 %`` is discouraged. Use `2.5 %` instead.
## Warning: Use of `df$`97.5 %`` is discouraged. Use `97.5 %` instead.
save plots
pdf(file = "analysisOutputs/plots/vr_ci_plot.pdf", width = 7, height = 5)
vr.ci.plot
## Warning: Use of `df$`2.5 %`` is discouraged. Use `2.5 %` instead.
## Warning: Use of `df$`97.5 %`` is discouraged. Use `97.5 %` instead.
dev.off()
## quartz_off_screen
## 2
pdf(file = "analysisOutputs/plots/cond_ci_plot.pdf", width = 7, height = 5)
cond.ci.plot
## Warning: Use of `df$`2.5 %`` is discouraged. Use `2.5 %` instead.
## Warning: Use of `df$`97.5 %`` is discouraged. Use `97.5 %` instead.
dev.off()
## quartz_off_screen
## 2
# save the plots
tikz(file = "analysisOutputs/plots/check.tex", width = 7, height = 5)
grid::grid.draw(shift_legend(vr_eval_plot))
## Warning: Removed 304 rows containing non-finite values (stat_ydensity).
## Warning: Removed 304 rows containing missing values (geom_point).
dev.off()
## quartz_off_screen
## 2
data.desc <- subset(data, subset =( vr & time==1) ) %>%
ungroup()
scls <- c("excitement", "graphically pleasing", "pleasant", "realistic", "enjoyment")
forms <- paste("vr_eval",1:5, " ~ condition", sep = "")
vr_evals <- list()
for(i in 1:5){
paste("vr_eval", i, " ~ condition", sep = "")
formul <- as.formula(paste("vr_eval", i, " ~ condition", sep = ""))
vr_evals[["lm"]][[paste("vr_eval", i, ": ", scls[i])]] <- lm.temp <- lm(data.desc, formula = formul)
vr_evals[["comparisons"]][[scls[i]]] <- anova(lm.temp)
}
vr_evals
## $lm
## $lm$`vr_eval 1 : excitement`
##
## Call:
## lm(formula = formul, data = data.desc)
##
## Coefficients:
## (Intercept) conditiona conditionc
## 6.34783 -0.17391 0.08696
##
##
## $lm$`vr_eval 2 : graphically pleasing`
##
## Call:
## lm(formula = formul, data = data.desc)
##
## Coefficients:
## (Intercept) conditiona conditionc
## 5.9565 -1.0435 -0.6087
##
##
## $lm$`vr_eval 3 : pleasant`
##
## Call:
## lm(formula = formul, data = data.desc)
##
## Coefficients:
## (Intercept) conditiona conditionc
## 6.0000 -0.1739 0.1739
##
##
## $lm$`vr_eval 4 : realistic`
##
## Call:
## lm(formula = formul, data = data.desc)
##
## Coefficients:
## (Intercept) conditiona conditionc
## 5.4348 -1.0000 -0.3478
##
##
## $lm$`vr_eval 5 : enjoyment`
##
## Call:
## lm(formula = formul, data = data.desc)
##
## Coefficients:
## (Intercept) conditiona conditionc
## 6.30435 -0.08696 0.17391
##
##
##
## $comparisons
## $comparisons$excitement
## Analysis of Variance Table
##
## Response: vr_eval1
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 0.812 0.4058 0.6667 0.5168
## Residuals 66 40.174 0.6087
##
## $comparisons$`graphically pleasing`
## Analysis of Variance Table
##
## Response: vr_eval2
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 12.638 6.3188 3.208 0.04682 *
## Residuals 66 130.000 1.9697
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $comparisons$pleasant
## Analysis of Variance Table
##
## Response: vr_eval3
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 1.391 0.69565 0.4307 0.6519
## Residuals 66 106.609 1.61528
##
## $comparisons$realistic
## Analysis of Variance Table
##
## Response: vr_eval4
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 11.855 5.9275 4.49 0.01485 *
## Residuals 66 87.130 1.3202
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $comparisons$enjoyment
## Analysis of Variance Table
##
## Response: vr_eval5
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 0.812 0.40580 0.552 0.5784
## Residuals 66 48.522 0.73518
summary(vr_evals$lm$`vr_eval 2 : graphically pleasing`)
##
## Call:
## lm(formula = formul, data = data.desc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3478 -0.9565 0.0870 1.0435 2.0870
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.9565 0.2926 20.354 <2e-16 ***
## conditiona -1.0435 0.4139 -2.521 0.0141 *
## conditionc -0.6087 0.4139 -1.471 0.1461
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.403 on 66 degrees of freedom
## Multiple R-squared: 0.0886, Adjusted R-squared: 0.06098
## F-statistic: 3.208 on 2 and 66 DF, p-value: 0.04682
summary(vr_evals$lm$`vr_eval 4 : realistic`)
##
## Call:
## lm(formula = formul, data = data.desc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0870 -0.4348 0.5652 0.5652 1.9130
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4348 0.2396 22.685 < 2e-16 ***
## conditiona -1.0000 0.3388 -2.951 0.00438 **
## conditionc -0.3478 0.3388 -1.027 0.30836
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.149 on 66 degrees of freedom
## Multiple R-squared: 0.1198, Adjusted R-squared: 0.09309
## F-statistic: 4.49 on 2 and 66 DF, p-value: 0.01485
Here we see that only the realism was rated clearly differently between the conditions, although excitement and graphical pleasantness are borderline.
vars <- c("sod", "ipq")
forms <- paste(vars, " ~ condition", sep = "")
sod_ipq_list <- list()
for(i in 1:2){
formul <- as.formula(forms[i])
sod_ipq_list[["lm"]][[vars[i]]] <- lm.temp <- lm(data.desc, formula = formul)
sod_ipq_list[["comparisons"]][[scls[i]]] <- anova(lm.temp)
}
sod_ipq_list
## $lm
## $lm$sod
##
## Call:
## lm(formula = formul, data = data.desc)
##
## Coefficients:
## (Intercept) conditiona conditionc
## 4.888889 0.280193 -0.009662
##
##
## $lm$ipq
##
## Call:
## lm(formula = formul, data = data.desc)
##
## Coefficients:
## (Intercept) conditiona conditionc
## 4.7112 -0.2554 0.1324
##
##
##
## $comparisons
## $comparisons$excitement
## Analysis of Variance Table
##
## Response: sod
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 1.247 0.62337 0.5808 0.5623
## Residuals 66 70.834 1.07324
##
## $comparisons$`graphically pleasing`
## Analysis of Variance Table
##
## Response: ipq
## Df Sum Sq Mean Sq F value Pr(>F)
## condition 2 1.635 0.81748 1.632 0.2038
## Residuals 62 31.057 0.50092
Now comparing vr to non-vr
vars <- c("sod", "ipq")
forms <- paste(vars, " ~ vr", sep = "")
sod_ipq_vr_list <- list()
for(i in 1:2){
formul <- as.formula(forms[i])
sod_ipq_vr_list[["lm"]][[vars[i]]] <- lm.temp <- lm(data, formula = formul, subset = time == 1)
sod_ipq_vr_list[["comparisons"]][[scls[i]]] <- anova(lm.temp)
}
sod_ipq_vr_list
## $lm
## $lm$sod
##
## Call:
## lm(formula = formul, data = data, subset = time == 1)
##
## Coefficients:
## (Intercept) vrTRUE
## 4.6250 0.3541
##
##
## $lm$ipq
##
## Call:
## lm(formula = formul, data = data, subset = time == 1)
##
## Coefficients:
## (Intercept) vrTRUE
## 3.586 1.085
##
##
##
## $comparisons
## $comparisons$excitement
## Analysis of Variance Table
##
## Response: sod
## Df Sum Sq Mean Sq F value Pr(>F)
## vr 1 2.232 2.2323 1.886 0.173
## Residuals 91 107.706 1.1836
##
## $comparisons$`graphically pleasing`
## Analysis of Variance Table
##
## Response: ipq
## Df Sum Sq Mean Sq F value Pr(>F)
## vr 1 20.639 20.639 37.323 2.739e-08 ***
## Residuals 87 48.110 0.553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
data.desc.all <- subset(data, subset = time==1) %>%
ungroup()
scls <- c("excitement", "graphically pleasing", "pleasant", "realistic", "enjoyment")
forms <- paste("vr_eval",1:5, " ~ condition", sep = "")
vr_evals_all <- list("lm" = list(), "comparisons" = list())
for(i in 1:5){
#i <- 1
paste("vr_eval", i, " ~ condition", sep = "")
formul <- as.formula(paste("vr_eval", i, " ~ vr + (1|condition)", sep = ""))
vr_evals_all[["lm"]][[scls[i]]] <- mem.temp <- lmer(data.desc.all, formula = formul)
vr_evals_all[["comparisons"]][[scls[i]]] <- anova(mem.temp)
}
## boundary (singular) fit: see ?isSingular
summary(lm(data.desc.all, formula = ipq ~ vr))
##
## Call:
## lm(formula = ipq ~ vr, data = data.desc.all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.65774 -0.52857 -0.02857 0.54286 1.62798
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.5863 0.1518 23.626 < 2e-16 ***
## vrTRUE 1.0851 0.1776 6.109 2.74e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7436 on 87 degrees of freedom
## (53 observations deleted due to missingness)
## Multiple R-squared: 0.3002, Adjusted R-squared: 0.2922
## F-statistic: 37.32 on 1 and 87 DF, p-value: 2.739e-08
summary(lm(data.desc.all, formula = sod ~ vr))
##
## Call:
## lm(formula = sod ~ vr, data = data.desc.all)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.73611 -0.53462 -0.09018 0.70833 2.37500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.6250 0.2221 20.827 <2e-16 ***
## vrTRUE 0.3541 0.2578 1.373 0.173
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.088 on 91 degrees of freedom
## (49 observations deleted due to missingness)
## Multiple R-squared: 0.0203, Adjusted R-squared: 0.009539
## F-statistic: 1.886 on 1 and 91 DF, p-value: 0.173
Looks like VR has an impact on enjoyment and on excitement. It seems to be a positive impact! VR also leads to larger presence, but not to larger suspension of disbelief than video.